      subroutine readinput

      use param
      use allocs

#include <global.h>
      include 'geom.h'
      include 'geomf.h'
      include 'xsec.h'
      include 'cmfd.h'
      include 'cmfdp1.h'

      localfn='sample_2g.inp'
!      localfn='c5g7.inp'
!      localfn='moc_g7.inp'
      open (unit=71,file=localfn,status="old")

      read(71,*) nreg,ncomp,ng

      call dmalloc(rgz,nreg)
      call dmalloc(acom,nreg)
      call dmalloc(sn,nreg)
      call dmalloc(pn,nreg)
      call dmalloc(diff,ncomp,ng)
      call dmalloc(sigr,ncomp,ng)
      call dmalloc(signf,ncomp,ng)
      call dmalloc(chi,ncomp,ng)
      call dmalloc(sigt,ncomp,ng)
      call dmalloc(rsigt,ncomp,ng)

      allocate (sigsm(ncomp,ng))
      
      do im=1,ng
        do icomp=1,ncomp
          allocate(sigsm(icomp,im)%from(ng))
        enddo
      enddo

      read(71,*) rgz(1:nreg)
      read(71,*) acom(1:nreg)
      read(71,*) sn(1:nreg)
      read(71,*) pn(1:nreg)
      read(71,*) opt
      read(71,*) bcb,bcu
      read(71,*) kconv,fconv
      read(71,*) nout

      tnode=sum(sn)    
      fnode=sum(pn)

      call dmalloc(fn,tnode+1)
      pe=0
      do i=1,nreg
        sb=1
        se=se+sn(i)
        do j=sb,se
          fn(j)=pn(i)/sn(i)   
        enddo
      enddo
      fn(tnode+1)=fn(tnode)
 
      call allocpdm

!     separation
      k=1
      do j=1,nreg
        do i=k,k+sn(j)-1
          hz(i)=rgz(j)/sn(j)
          rhz(i)=1/hz(i)
          comp(i)=acom(j)
        enddo
        k=i
      enddo

      k=1
      do j=1,nreg
        do i=k,k+pn(j)-1
          hzf(i)=rgz(j)/pn(j)
          rhzf(i)=1/hzf(i)
          compf(i)=acom(j)
        enddo
        k=i
      enddo

!     read xsec
      im=1
      icomp=1
      do while (icomp.lt.ncomp)
        read(71,*) icomp
        im=1
        do while (im.lt.ng)
          read(71,*) im,diff(icomp,im),sigr(icomp,im),signf(icomp,im) &
     ,              chi(icomp,im)
        enddo
        im=1
        do while (im.lt.ng)
          read(71,*) im,(sigsm(icomp,im)%from(ms),ms=1,ng)
        enddo
      enddo

!     define total xsec
      do im=1,ng
        do icomp=1,ncomp
          sigt(icomp,im)=sigr(icomp,im)+sigsm(icomp,im)%from(im)   ! removal+self scattering 
          rsigt(icomp,im)=1/sigt(icomp,im)         
          sigsm(icomp,im)%from(im)=0        ! set 0 for self-scattering  
        enddo
      enddo

      end subroutine
